clear;clc;
A=[0.98 0.1 0;0.02 0.7 0.05;0 0.2 0.95];
T=1000000;
states=zeros(1,T);
[b,v]=eig(A);

%(a)
%State C1 is represented by 1,state C2 by 2 and state O by 3,respectively
states(1)=2;
for t=1:T-1
   states(t+1)=myfun(states(t),A);
end
figure;
hist(states,6);
title('the distribution of states');
xlim([0 4]);
set(gca,'XTick',[1 2 3]);

%(b)
a1=length(find(states==1));
a2=length(find(states==2));
a3=length(find(states==3));
A1=1;A2=a2/a1;A3=a3/a1;
B1=1;B2=b(2,2)/b(1,2);B3=b(3,2)/b(1,2);
figure;plot([1,2,3],[A1 A2 A3],'o',[1 2 3],[B1 B2 B3],'+');
axis([0 4 0 1.5]);
legend('eigenvector','states','Location','best');
title('comparasion between states and dominant eigenvector with relative value');

%(c)
%State "open" is represented by -0.5,state "close" by 0.5,respectively
rstates=zeros(1,T);
for i=1:T
    if states(i)<3
        rstates(i)=-0.5;
    else
        rstates(i)=0.5;
    end;
end;
figure;
hist(rstates,5);
title('the distribution of rstates');
xlim([-1 1]);
set(gca,'XTick',[-0.5 0.5]);

%(d)
open=zeros(1,T);
closed=zeros(1,T);
 m=1;
 jo=1;
 jc=1;
for k=2:T
    if (rstates(k)==rstates(k-1))
        m=m+1;
    else  
       if (rstates(k-1) == 0.5)
           closed(jc)=m;
           jc=jc+1;
       else
           open(jo)=m;
           jo=jo+1;
       end;
       m=1;
    end
end;
open((open==0))=[];
closed((closed==0))=[];

figure;hist(open,0:1:1000);
xlabel('residence time');
ylabel('number of open');
axis([0 1000 0 5000]);
figure;hist(closed,0:1:1000);
xlabel('residence time');
ylabel('number of closed');
axis([0 1000 0 1000]);

%(e)
% the closed state is not fit well by an exponential distribution
figure;hist(closed,0:1:1000);
xlabel('residence time');
ylabel('number of closed');
axis([0 1000 0 1000]);
hold on;
% the exponential distribution
bTotal=b(1,2)+b(2,2)+b(3,2);
x=1:1:1000;
y=10000./(x.^bTotal);
plot(x,y,'r');

% the open state seems to fit well by an exponential distribution
figure;hist(open,0:1:1000);
xlabel('residence time');
ylabel('number of open');
axis([0 1000 0 5000]);
hold on;
% the exponential distribution
bTotal=b(1,2)+b(2,2)+b(3,2);
x=1:1:1000;
y=10000./(x.^bTotal);
plot(x,y,'r');